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Abstract 

We study the problem of computing outer bounds for the region of 
steady states of biochemical reaction networks modelled by ordinary dif- 
ferential equations, with respect to parameters that are allowed to vary 
within a predefined region. Using a relaxed version of the correspond- 
I ing feasibility problem and its Lagrangian dual, we show how to compute 

O"* certificates for regions in state space not containing any steady states. 

Based on these results, we develop an algorithm to compute outer bounds 
I for the region of all feasible steady states. We apply our algorithm to 

the sensitivity analysis of a Goldbeter-Koshland enzymatic cycle, which 
is a frequent motif in reaction networks for regulation of metabolism and 
^\ signal transduction. Copyright © 2008 IF AC. 

in 

^j.' 1 Introduction 

O 

Q\ A basic question in the analysis of biochemical reaction networks is how steady 

<^> state concentrations change with parameters. Metabolic Control Analysis (MCA) 



X 



is a classical tool to answer this question Kacser et al. 1995 , where the analysis 
is based on a linear approximation of the system's equations around the steady 



state. Due to the linear approximation, results from MCA are only valid if pa- 
rameter variations are small. However, in natural biochemical reaction networks, 
one usually faces large parameter variations: in genetic engineering, common 
techniques like gene knock-outs or knock-downs, overexpression or binding site 
mutations typically give rise to large parameter variations. 

It follows that there is a need to compute changes in steady state values 
which are due to large parameter variations. One approach to broaden the 
validity of results from MCA to larger parameter variations is to include higher 
order approximations at the nominal point [Streif et al. 2007 . Although such 



an approach may extend the validity of the approximation, it still gives results 
which are in general only locally valid. 
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We will thus rather take a different route and study the problem from the 
perspective of computing the set of all steady states for given ranges in which 
parameter values may vary. In contrast to classical, local sensitivity analysis, 
such an approach allows to directly evaluate the range that steady state con- 
centrations can take for given parameter ranges. The drawback is that it is not 
directly possible to assess the influence of individual parameters on the steady 
state. However, by repeating the computation for different parameter ranges, 
also this information may be obtained. 

Computing the set of steady states analytically is only possible in very rare 
cases. Even if an analytical solution for the steady state is known, computing 
the corresponding set for all possible parameter values may be difficult. Due 
to this difficulty, non-deterministic approaches are frequently used to solve this 
problem. A common tool for this kind of analysis are Monte Carlo methods 

which are routinely applied in the analysis of un- 



Robert and Casella 2004 



certain biochemical reaction networks Alves and Savageau 2000 Feng et al 



2004 . However, Monte Carlo methods do not give reliable results in the sense 
that it is possible to miss important solutions, which is particularly problem- 
atic for highly nonlinear dependencies of the steady state on parameters. Also, 
Monte Carlo approaches to the problem at hand typically require that all of the 
possibly multiple steady states for specific parameter values can be computed 
explicitly, which is often a difficult task in itself. 

Continuation methods that track the changes in steady state values upon 



parameter variations are an efficient computational tool for this problem Richtcr 
and DeCarlo 1983 Kuznetsov 1995 , but are restricted to low-dimensional 
parameter variations and are thus in general unsuitable for exploring higher- 
dimensional parameter spaces. 

Global optimization methods employing branch and bound techniques or in- 
terval arithmetics would in principle be suited to compute steady state regions 
Maranas and Floudas 1995 Ncumaier 1990 . However, it seems that the cor- 



responding computational cost has obstructed their application to the analysis 
of biochemical reaction networks so far. 

In this paper, we propose a new approach to obtain reliable bounds on steady 
state values under uncertain parameters in a computationally efficient way. The 
paper is structured as follows. In Section [3j we study the problem of finding 
certificates that a given set in state space does not contain a steady state for any 
parameters in a given set in parameter space. In Section [4j we use the results 
obtained in Section [3] to develop an algorithm that computes outer bounding 
regions of steady state values for a given set in which parameters vary. The 
application of the proposed analysis method is shown for two example reaction 
networks in Section [H 



Mathematical notation 

The space of real symmetric nxn matrices is denoted as S n . The order operator 
with respect to the positive orthant in R mxn is denoted as "<", i.e. 0<le 
K rax ™ <^> < Xij for i — l,...,m, j = l,...,n. The order operator with 
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respect to the cone of positive semidcfinitc (PSD) matrices in S n is denoted as 
V , i.e. X £ S n & X is PSD. The trace of a quadratic matrix X £ R" xrl 
is denoted as tlX. 



2 Problem statement and basic idea 

We consider biochemical reaction networks that are modelled by ordinary dif- 
ferential equations. This modelling framework is quite general and covers most 
metabolic networks as well as many signal transduction pathways, if spatial 
effects can be neglected. Mathematically, such models are commonly written as 

x = Sv(x,p), (1) 

where x £ K n is the concentration vector, S £ l" x ™ is the stoichiometric 
matrix, p £ K m is the vector of parameter values and v(x,p) £ K m is the vector 



of reaction fluxes Klipp et al. 2005 . Throughout this paper, we assume that 



fluxes are modelled using the law of mass action, where v takes the form 

n 

Vj{x,p) =Pj Y[x a k 3 \ (2) 

k=l 

for j = 1, . . . , m. The constants <Tjk are integers representing the stoichiometric 
coefficient of the species k taking part in the j'-th reacting complex. In the 
case of mass action kinetics, the dimensions of the parameter vector and the 
flux vector are in general the same. Note that our results can be extended 
to rational functions describing the fluxes, such as used for Michaelis-Menten 
kinetics, in a straightforward way. 

The problem under consideration can be formulated as follows. Given a set 
V C M. m in parameter space, compute a set X s C R" that contains all steady 
states of the system (JTJ for parameter values taken from V . Ideally, the set X s 
should be as small as possible, such that for all x s £ X s , there is a parameter 
vector p £ V with Sv(x s ,p) — 0. Then, 

X s = {x£R rn \3 P £V : Sv(x,p) = 0}. (3) 

However, for the case m > 1, when continuation methods are not suitable, there 
are at present no general methods to compute X s efficiently and reliably. 

We present a method to address this problem that works for arbitrarily 
large state and parameter spaces, does not need to compute steady state values 
explicitly and is computationally efficient. The method is able to compute 
reliable, though conservative outer bounds on the set X s of all steady states. 

In order to search for sets of steady states for a given parameter set V, we 
need means to test whether a candidate solution X s obtained in such a search is 
actually valid or not. Such a test is readily formulated as a feasibility problem. 
Moreover, we will see that the Lagrangian dual for this feasibility problem allows 
to certify given regions in state space as not containing a steady state for any 
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parameter value from the set V . We then develop an algorithm that uses this 
information to construct outer bounds on the region X s of all steady states. 

In this paper, we consider only hyperrectangles for the sets X s and V in 
state and parameter space. An extension to more general convex polytopes is 
in principle easy from the theoretical perspective, but it requires a much more 
elaborate implementation on the practical side. 



3 Feasibility of steady state regions 

3.1 Feasibility problem and semidefinite relaxation 

The problem of testing whether a given hyperrectangle X s in state space con- 
tains steady states of the system ([!]), for some parameter values in a given 
hyperrectangle V in parameter space, can be formulated as the following feasi- 
bility problem: 




, m 
,n. 



(4) 



The same problem appears in the context of parameter identification in a 
recent paper by Kuepfer et al. 2007| . They developped a method that uses an 
infeasibility certificate for the problem Q to exclude regions in parameter space 
from the identification procedure, given a set of steady state measurements. In 
this section, we take their approach to find an infeasibility certificate for prob- 
lem Q, but give more details about the underlying mathematical techniques. 

Relaxing the feasiblity problem Q to a semidefinite program Vandenberghe 
and Boyd 1996 ensures computational efficiency. The applied relaxation is 
based on a quadratic representation of a multivariate polynomial of arbitrary 
degree Parrilo 2003 . In the first step, we construct a vector £ containing 
monomials that occur in the reaction flux vector v(x,p). In the special case 
where no single reaction has more than two reagents, a starting point for the 
construction of £ is 

£ (1) Pi j • • • 3 Pra 3 <B\ i • • • ) X n j Pl%l 3 ■ ■ ■ ) Pj %i 3 ■ ■ ■ 3 Pm-En) : 

which can usually be reduced by eliminating components that are not required 
to represent the reaction fluxes. We define k such that £ € R fe . Note that 
this approach is not limited to second order reaction networks. In more general 
cases, one has to extend the vector £ by monomials that are products of several 
state variables. 

Using the vector £, the elements of the flux vector v(x,p) can be expressed 

as 

Vj(x,p) = £ T Vj€, i = l,...,m, (5) 
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where Vj G S is a constant symmetric matrix. The choice of Vj is generally 
not unique, as an expression of the form pjXiXk can be decomposed as either 
(j>jXi)(xk) or (pjXk)(xi). This fact may be used to introduce additional equality 
constraints in the relaxed problem ([8]), but we will neglect this for simplicity of 
notation. 

Using ([5]), the system ([I]) can be written as 



.,n, 



(6) 



where Qi = X^fLi G S k are constant symmetric matrices. 

The original feasibility problem Q is thus equivalent to the problem 



find 
s.t. 





B£ > 
Ci = l, 



(7) 



where the matrix £? G lR ( - 2fc 2 - )xfc is constructed to cover the inequality con- 
straints in Q, e.g. the constraint P\ >m i n < pi < Pi,mox is represented as 



-Pl.mtn 1 ■ ■ ■ 
Pl,max "I ... 



Corresponding constraints for higher order monomials in £ arc obtained easily 
as Pj, m i n %i,min < Pj%i < Pj,max%i,max and have to be included in the matrix B. 

A relaxation to a semidefinite program is found by setting X = ££ T . The re- 
sulting non-convex constraint rankX = 1 is omitted in the relaxation. Instead, 
several consequences of how X is defined, namely Xu = 1 and X )p 0, are used 
as convex constraints. The relaxed version of the original feasibility problem Q 
is thus obtained as 



(RP) 



find 
s.t. 



x e s k 

tr(QiX) = 
tr( ei eJX) = 1 
BX ei > 
BXB T > 

x^o, 



i = l,. 



,n 



(8) 



where e x = (1,0, .. .,0) T € M fc . 

The basic relationship between the original problem (jij and the relaxed 
problem |8| is that if the original problem is feasible, then the relaxed problem 
is also feasible. Thus, the relaxation allows to certify a region in state space 
as infeasible for steady states, as we will see when going to the Lagrange dual 
problem. 
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3.2 Infeasibility certificates from the dual problem 

The Lagrange dual problem can be used to certify infeasibility of the primal 
problem (|8]). First, the Lagrangian function L is constructed for the primal 
problem. We obtain 

L(X, Ai, A 2 , A 3 , v) = -\jBX ei - tr(Aj BXB T ) 

n 

- tr(XjX) + ^2^ti(QiX) + v n+1 (tr( ei eJX) - 1), 

i=l 

where Ai G R 2k - 2 , A 2 G S 2k - 2 , A 3 G S k and v G M" +1 . Using the cyclic 
property of the trace operator, i.e. tr(ABC) = tr(BCA) = tr(CAB), we rewrite 

tr(A 2 r BXB T ) = tr(B T X^BX) 

and 

XjBX ei = ti(e!-±-BX) +tr(ef'yS T X) 

= tr(( ei ^B + e ^B T )X). 

The second reformulation has also the advantage of providing a symmetric mul- 
tiplier for X, which is more efficient from the computational side. 
Based on the Lagrangian L, the dual problem is obtained as 

max inf L(X, Ai, A 2 , A 3 , u) 

xes k 

s.t. Ai > 0, A 2 > 0, A 3 0, 
which is equivalent to 

max f n +i 

s.t. B T X 2 B + e 1 XjB + B T X 1 eJ 



(D) 



+X 3 + ^2 "iQi + v n+ ieiej = 

i=l 

Ai > 0, A 2 > 0, A 3 ^ 0. 



(9) 



It is a standard procedure in convex optimisation to use the dual problem 
in order to find a certificate that guarantees infeasibility of the primal problem 



Boyd and Vandenberghe 2004 . For the problem at hand, this principle is 



formulated in the following theorem. 

Theorem 1. If the dual problem ^ has a feasible solution where v n+ \ > 0, 
then the primal problem Q is infeasible. 

Proof. Note that the constraints of the dual problem ^ are homogenous in the 
free variables: if (A' 1; A 2 , A 3 , v') is feasible, then also (aA^, aX' 2 , aX' 3 , av') with 
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any a > is feasible. In particular, choosing all free variables to be zero is 
always a feasible solution of the dual problem ^ . 

Let d* be the optimal value of the dual problem ([9| . By the previous argu- 
ment, it is clear that either d* = or d* = oo. Under the assumption made in 
the theorem, we have d* = oo. 

To the primal feasibility problem ^ , we can associate a minimization prob- 
lem with zero objective function and the same constraints as in pi). Let p* be 
the optimal value of this minimization problem. We have p* = 0, if the primal 
problem ([8| is feasible, and p* — oo otherwise. Weak duality of semidefinite 



programs Vandenberghe and Boyd 1996 assures that d* < p* . In particular, 



d* = oo implies p* = oo, and the primal problem ([8| as well as the original 
feasibility problem Q are both infeasible. □ □ 

Theorem [I] sets the basis for our further considerations. 



4 Bounding feasible steady states 

In this section, we present an approach to find bounds on the steady state region 
X s , based on the results obtained in the previous section. As basic additional 
requirement, we assume that some upper and lower bounds on steady states are 
already known previously by other means. Let these bounds be given by 

•Ei, lower ^ Xi — -^i,upper: % — 1, . . . , 77-. (10) 

In biochemical reaction networks, such bounds can often be obtained from mass 
conservation relations, as done for the examples in Section [5j Also, it is often 
possible to show positive invariance of a sufficiently large compact set in state 
space for the system (JT|). These bounds may be very loose though, and the main 
objective of our method is to tighten them as far as possible. 

To this end, we use a bisection algorithm that finds the maximum ranges 
[x j dower, Xj,min] and [xj^max , %j, upper) for which infcasibility can be proven via 
Theorem [l] The algorithm iterates over j — 1, . . . ,n, while the steady state 
values Xi for i =/= j are assumed to be located within the interval given by 
inequality ( |To| ). 

We give the bisection algorithm in pseudocode for computing the lower 
bound xi^min. The computation of the upper bound xi^ max works in essen- 
tially the same way, with some obvious modifications. 

Algorithm 1 (Lower bound maximization by bisection). 

Up_gUeSS <- X X ,upper 
l0_gUeSS <- Xx^ower 
nGX"t_X]. ^ — X\^ upper 

while (up_guess - lo_guess) > tolerance 
use constraint X\j ower < X\ < next_Xi 
solve semidefinite program (D) 

if d* = oo 
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lo_guess <- next_a;i 

increase next_.xi by i(up_guess - next_a;i) 
else 

up_guess <- next_a;i 

decrease next_xi by i(next_xi ~~ lo_guess) 
endif 
endwhile 

%l,min <~ lo_guess 

Due to the availability of efficient solvers for semidefinite programs and the 
use of bisection to maximize the interval that is certified as infeasible, Algo- 
rithm [T] can run considerably fast on standard desktop computers, as we will 
see in the examples discussed in the following section. 

In our analysis method, Algorithm [T] is run for all state variables, and as both 
maximization of the lower bound and minimization of the upper bound of the 
steady state values. Its output is a hyperrectangle in state space containing all 
steady states for the assumed parameter ranges. This is a relevant information 
for the global sensitivity analysis of a biochemical reaction network, as it allows 
to discriminate concentration values that arc highly affected by the assumed 
parameter variations from others that are less affected. Moreover, by repeating 
the computation for different parameter ranges, it is also possible to assess the 
influence of individual parameters on steady state concentrations, which is closer 
related to classical, local sensitivity analysis. 



5 Examples 

5.1 A simple conversion reaction 

As first example, we consider a simple conversion reaction where the region of 
steady states for a given parameter box can be computed analytically. Consider 
the reaction network 

fei 

A <=i B. 

Denote the concentrations of A and B as a and 6, respectively. There is a 
conservation relation a(t) + b(t) — clq, so the system can be modelled by one 
differential equation 

a = k2(ao — a) — k\a. (11) 

Furthermore, there is a unique steady state a s for all parameter values, given 
by 

_ k 2 ao 
ki + k 2 

From the conservation relation, we have the loose bound < a s < ao which is 
valid for all parameter values. Assume now that ao = 1 is fixed, and let the 
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other parameters vary in a box k\,k 2 € [k m i n , k max ]. Then, the steady state 
varies in the interval 

h h 

"'mm ^max 
as h 4- h ■ ' h 4- h 

In the specific case where k m i n = 1 and k max = 2, the steady state interval is 
a s G [gjf]- Our algorithm is able to compute numerically exact bounds in these 
cases. For a numerical precision of 10~ 6 , computation time is a few seconds on 
a standard desktop computer. 



5.2 An enzymatic cycle 



As a more complex example, where the steady state region for a given parameter 
box cannot be computed analytically, we consider an enzymatic cycle. These 
cycles appear very frequently in cellular reaction networks, in particular in the 
form of phosphorylation/dephosphorylation cycles |Shacter et al. 1984 . An 



enzymatic cycle as encountered in covalcnt modification of proteins Goldbeter 



and Koshland 1981 is typically described by the reaction network 

fei 



E + A 



Ci 



P 



Ci 
A* 

Co 



E + A* 



(12) 



Co 



fer, 



P + A. 



There are three conservation relations 

[A] + [A*] + [d] + [C 2 ] = A Q 
[E] + [Ci] = E 
[P] + [Ch] = Po- 

Denoting a = [A*], C\ = [Ci] and c-i — [Co] and using the law of mass action, 
the reaction flux vector is given by 

fki(A -a-ci- c 2 )(E - ci)\ 
k 2 ci 

fc 3 Ci 

k 4 (P Q - c 2 )a 

k 5 C2 

V hc 2 J 

Due to the conservation relations, we only need to use three differential equations 
in the model, which is given by 



dt 



(a, ci, c 2 ) T = 




(13) 
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For the sensitivity analysis, the parameters k\ and k± as well as the total con- 
centrations Aq, E and P are assumed to be fixed at fei = 10 5 , fc 4 = 5 • 10 4 , 
Ao = 1 and Eq = Po = 0.01. The other parameters are assumed to be variable 
parameters, with variations around their nominal values fc 2j „ m = k 5 nom = 1 
and k 3nom ^6,nom 10 ■ 

From the conservation relations and invariance of the positive orthant we 
have the steady state bounds 

< a < Aq, < ci < E , < c 2 < P , 

which are valid for any parameter values. 

We have applied the proposed analysis method to find tighter bounds on 
possible steady state values, comparing three different regions in which param- 
eters of the enzymatic cycle are allowed to vary. The three different regions are 
given by Vi, V 2 and V 3 , where Vi, V 2 , V 3 C M 4 and 

• (fc 2 , fc 3 , fc 5 , fc 6 ) e Vi <^ 0.98ki, nom < h < 1.02fc linom , corresponding to 
parameter variations of up to 2%, 

• (fc 2 , fc 3 , fc 5 , fc 6 ) £ P 2 0.9fc^„ om < fcj < l.lfcj )nom , corresponding to 
parameter variations of up to 10%, and 

• (fc 2 ,fc 3 ,fc 5 ,fc 6 ) G P 3 <^ 0.5fc ijnom < ki < 2ki. nom , corresponding to up to 
2-fold parameter variations, 

with i = 2,3,5,6 in all three cases. 

The dual problem (D) has been constructed by using 

£ T = (1, k 2 , h, k 5 , k 6 , a, a, c 2 ), (14) 

and deriving appropriate matrices Qi, B 7 for the steady state equations and 
the constraints, respectively. Algorithm [T| was then used to compute bounds 
on the steady state concentrations. We compare these results to an estimate 
for the region of steady state concentrations obtained by Monte-Carlo tests. 
The results are shown in Figure [l] The average computation time to obtain 
the feasible intervals for all three state variables and one parameter region was 
about 25 seconds. The Monte-Carlo tests done to produce the figures took 
consistently about 20 % more computation time, where 1000 parameter points 
were used for each test. However, for a reliable evaluation by Monte-Carlo 
methods, much more points should be used, which would increase computation 
time significantly. 

As can be seen from the figure, our approach is able to find tight intervals for 
the steady state values of the individual concentrations. However, the results 
also highlight the limitations of using hyperrectanglcs if the steady state values 
are highly correlated. 

Our analysis also yields a biochemical interpretation, related to the property 
of ultrasensitivity. The concept of ultrasensitivity is quite important for bio- 
chemical reaction networks, in particular for those that constitute cellular signal 
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0.2 0.4 0.6 0.8 
Steady state value of [A'] 



02 0.4 0.6 0.8 

Steady state value of [A*] 



(a) Parameter uncertainty region 
variation) 



Pi (2% (b) Parameter uncertainty region Vi (10% 
variation) 




0.2 0.4 0.6 0.8 

Steady state value of [A*] 



(c) Parameter uncertainty region Vz (2-fold 
variation) 



Figure 1: Feasible steady states for the enzymatic cycle with three different 
parameter regions, comparison of reliable bounds obtained with Algorithm [T] 
and Monte-Carlo estimates. Light gray regions have been certified infeasible 
by Algorithm [T] Black dots are steady state values obtained from Monte-Carlo 
tests. Dark gray regions are known to be infeasible from conservation relations. 
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p 



Figure 2: Illustration of ultrasensitivity. Response of an output variable a; to a 
control variable p. The response is ultrasensitive, with high sensitivity around 
the nominal value po and considerable less sensitivity for other values. 



transduction pathways Levine et al. 2007 . Shortly, ultrasensitivity means that 



a small variation in a control variable has a relatively large effect on an out- 
put variable, whereas for increasing variations in the control variable, the range 
of the output variable will be considerably less increasing (see also Figure [2]) . 
Thus, ultrasensitivity is an inherently non-linear and non-local property. For 
the enzymatic cycle, a variation of only 2% in parameters already allows the 
steady state value of [A*] to vary over almost half of the interval given from 
the conservation relation, and with an allowable parameter variation of 10% 
the steady state value of [A*] can span nearly the whole interval. This is a 
clear indication of the ultrasensitivity which is typical for the enzymatic cycle 
|Goldbeter and Koshlan"d| |1981| . 

In addition, our results show that the steady state value of [Ci], the con- 
centration of the intermediate enzyme-substrate complex, is not ultrasensitive, 
because its value spans a large interval only for large parameter variations. 
Similar results hold for [C2]. 



6 Conclusions 

We have studied the problem of computing the region of all steady states of 
biochemical reaction networks, provided that parameters are allowed to vary 
within a known region. This is an important problem in sensitivity analysis of 
reaction networks. Our approach is based on formulating a feasibility problem 
to check whether a candidate region in state space actually contains steady 
states. This feasibility problem is relaxed to a semidefinite program, and its 
Lagrangian dual provides certificates of infeasibility of a candidate region in 
state space. These certificates can be used to efficiently minimize the estimate 
of the known feasible region in state space by a bisection algorithm. 

We have applied our sensitivity analysis to two simple example networks. For 
the first example, our algorithm is able to compute numerically exact bounds, 
which could be verified from the analytical solution. In the second example, 
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we compared the bounds obtained from our algorithm to steady state values 
obtained through Monte-Carlo tests. In this example, our approach was more 
efficient computationally than Monte-Carlo tests. Also, it gives guaranteed 
bounds on the steady state values, which cannot be achieved by randomized 
methods such as Monte-Carlo tests. Based on the premise that we are working 
with hyperrectangles only, the obtained bounds are fairly tight. The second 
example also shows that our approach is able to confirm ultrasensitivity of the 
Goldbeter-Koshland switch. 

In summary, our approach is a reliable and computationally efficient method 
to estimate the range of possible steady state variations due to multiple simulta- 
neous parameter variations in biochemical reaction networks, and thus provides 
a valuable tool for global sensitivity analysis. 
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